function[] = QC_d1d2(data,T,alpha1,alpha2,alpha3,t1,t2L,t2U,t3,stderr)
%%%% Computes CDF-based (quick&crude) estimates of the change-points d1, d2
%%%% and their standard error based on bootstrap.
%%%% INPUT ARGUMENTS
%%%% data = vector of all bid times
%%%% T = length of auctions (in seven day auctions, enter 7)
%%%% alpha1 = estimate of alpha1 (you can use QC_alpha.m to get an estimate)
%%%% alpha2 = estimate of alpha2 (you can use QC_alpha.m to get an estimate)
%%%% alpha3 = estimate of alpha3 (you can use QC_alpha.m to get an estimate)
%%%% t1 = a value that is likely to be in the first (espresso) phase [0, d1]
%%%% t2L, t2U = values that are likely to be in the second (macchiato) phase [d1, T-d2], such that t2L < t2U
%%%% t3 = a value that is likely to be in the third (ristretto) phase [T-d2, T]
%%%% stderr = flag for whether bootstrap stanadard deviation is computed.
%%%% Default is yes. Chose value other than 1 to supress this.
if nargin < 10,
    stderr = 1;
end

[d1 d2] = QC_estimate_d1d2(data,T,alpha1,alpha2,alpha3,t1,t2L,t2U,t3)

if stderr == 1,
    M=500; %%% Set number of bootstrap samples
    bt_data = bootrsp(data,M);

    for j=1:M, 
        if QC_estimate_d1d2(bt_data(:,j)',T,alpha1,alpha2,alpha3,t1,t2L,t2U,t3) < Inf,
            [d1(j) d2(j)]= QC_estimate_d1d2(bt_data(:,j)',T,alpha1,alpha2,alpha3,t1,t2L,t2U,t3); 
        end
    end
    std_d1 = std(d1)/sqrt(M)
    std_d2 = std(d2)/sqrt(M)
end

function[d1,d2] = QC_estimate_d1d2(data,T,alpha1,alpha2,alpha3,t1,t2L,t2U,t3)
d1 =  T - T* (  alpha1/alpha2 * Empirical_CDF(data,t1)/(Empirical_CDF(data,t2U)-Empirical_CDF(data,t2L)) *...
    ( (1-t2L/T)^alpha2-(1-t2U/T)^alpha2 )/ (1-(1-t1/T)^alpha1) )^(1/(alpha2-alpha1));

d2 = T* (  alpha3/alpha2 * (1-Empirical_CDF(data,t3))/(Empirical_CDF(data,t2U)-Empirical_CDF(data,t2L)) *...
    ( (1-t2L/T)^alpha2-(1-t2U/T)^alpha2 )/ (1-t3/T)^alpha3 )^(1/(alpha2-alpha3));
